Analysis of stimulation experiments.

library(Seurat)
## Loading required package: SeuratObject
## Loading required package: sp
## 
## Attaching package: 'SeuratObject'
## The following objects are masked from 'package:base':
## 
##     intersect, saveRDS
## Loading Seurat v5 beta version 
## To maintain compatibility with previous workflows, new Seurat objects will use the previous object structure by default
## To use new Seurat v5 assays: Please run: options(Seurat.object.assay.version = 'v5')
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(ggplot2)

Load integrated object. This object has been integrated with k.anchor = 20, because there was lots of batch effect between the 0hr tissue and the other two experiments (which merged quite well).

load("~/Dropbox/Zoe/scf_version/analysis/healthy_sc/seurat_objects/no_dropletQC/3391/integrated_3391_stimulation_pmaio_unstim_downsampled_cca_kanchor5_scClustViz.RData")

Set identities:

Idents(scSeurat) <- "integrated_snn_res.0.4"

Check out the UMAP of clusters, original identities, and refined SCINA labels.

DimPlot(scSeurat, label = TRUE)

DimPlot(scSeurat, group.by = "tissue_type", label = FALSE)

DimPlot(scSeurat, group.by = "scina_labels", label = TRUE) & NoLegend()

Look at tissue-type per cluster on a proportional level:

pt <- table(Idents(scSeurat), scSeurat$tissue_type)
pt <- as.data.frame(pt)
ggplot(pt, aes(x = Var1, y = Freq, fill = Var2)) +
  theme_bw(base_size = 15) +
  geom_col(position = "fill", width = 0.5) +
  xlab("Sample") +
  ylab("Proportion") +
  theme(legend.title = element_blank())

Look for T-cell signatures.

FeaturePlot(scSeurat, features = c("CD3E", "NKG7"))

FeaturePlot(scSeurat, features = c("LEF1", "XCL1;XCL2"))

FeaturePlot(scSeurat, features = c("TOX", "CD8A"))

Look for macrophage signatures:

FeaturePlot(scSeurat, features = c("LGALS1", "TYROBP"))

FeaturePlot(scSeurat, features = c("MARCO", "C1QC"))

FeaturePlot(scSeurat, features = c("LYZ-1", "S100A8"))

FeaturePlot(scSeurat, features = c("DEFA6;DEFA4;DEFA5", "CTSS"))

FeaturePlot(scSeurat, features = c("FLT3", "S100A12"))

Look for B and plasma cell signatures

FeaturePlot(scSeurat, features = c("CD79B", "JCHAIN"))

Cluster 12 is the only cluster that is obviously dominated by the unstimulated experiments. What cells are these? Looks like they are hepatocytes. A lot of the markers are periportal, but all hepatocytes are in this cluster.

# Isolate DE genes
clust11 <- sCVdata_list$res.0.4@DEvsRest$`12`[order(sCVdata_list$res.0.4@DEvsRest$`12`$Wstat,
                                                    decreasing = TRUE), ]
rownames(clust11)[1:50]
##  [1] "PAH"                                     
##  [2] "LOC114082496"                            
##  [3] "TDO2"                                    
##  [4] "MAT1A"                                   
##  [5] "GYS2"                                    
##  [6] "SLCO1B3;SLCO1B3-SLCO1B7;SLCO1B7;SLCO1B1" 
##  [7] "FBP1"                                    
##  [8] "CPS1"                                    
##  [9] "SLC2A2"                                  
## [10] "ALDOB"                                   
## [11] "MTSS1"                                   
## [12] "IGFBP1"                                  
## [13] "CYP3A4;CYP3A7-CYP3A51P;CYP3A7-2"         
## [14] "mikado.WCK01-AAH20201022-F8-rsc00143G668"
## [15] "PCK1"                                    
## [16] "SHMT1"                                   
## [17] "CYP3A4;CYP3A7-CYP3A51P;CYP3A7"           
## [18] "ARG1"                                    
## [19] "CYP3A4;CYP3A7-CYP3A51P;CYP3A7-3"         
## [20] "LOC107139912"                            
## [21] "ALB-1"                                   
## [22] "SLC38A4"                                 
## [23] "LOC107139914"                            
## [24] "SORBS2"                                  
## [25] "TAT"                                     
## [26] "C5-1"                                    
## [27] "KLKB1"                                   
## [28] "Slc22a19-1"                              
## [29] "G6PC1"                                   
## [30] "RBP4"                                    
## [31] "AGMO"                                    
## [32] "HSD17B11"                                
## [33] "CMBL"                                    
## [34] "SLC25A13"                                
## [35] "ZBTB16"                                  
## [36] "C8G"                                     
## [37] "PYGL"                                    
## [38] "HSD11B1"                                 
## [39] "mikado.WCK01-AAH20201022-F8-rsc00060G101"
## [40] "CES3"                                    
## [41] "MTARC2"                                  
## [42] "IDO2"                                    
## [43] "AQP9"                                    
## [44] "ACACB"                                   
## [45] "GOT1"                                    
## [46] "CFI"                                     
## [47] "NR1I2"                                   
## [48] "ABCB11"                                  
## [49] "SOX5"                                    
## [50] "BHMT2"
FeaturePlot(scSeurat, features = c("PAH","PCK1"))

FeaturePlot(scSeurat, features = c("CYP2E1","FETUB"))

FeaturePlot(scSeurat, features = c("ACACB","ELOVL6"))

FeaturePlot(scSeurat, features = c("POLR2D", "ALB-1"))

Let’s look at some of the key markers that Sonya sent: IL2, IL4, IL6, IL7, IL10, IL12A, TGFB

FeaturePlot(scSeurat, features = c("IL2", "IL4"))

FeaturePlot(scSeurat, features = c("IL6", "IL6-1"))

FeaturePlot(scSeurat, features = c("IL7", "IL10"))

FeaturePlot(scSeurat, features = c("IL12A", "TGFB1"))

It’s pretty hard to see which cytokines are expressed in which samples, so I’m going to try splitting up the samples.

DimPlot(scSeurat, split.by = "tissue_type", label = TRUE) & NoLegend()

Now let’s take a look at the features on this split plot.

FeaturePlot(scSeurat, feature = "IL2", split.by = "tissue_type") # Only in PMAIO cluster 2
## Warning: All cells have the same value (0) of "IL2"

FeaturePlot(scSeurat, feature = "IL4", split.by = "tissue_type") # Only in PMAIO cluster 2
## Warning: All cells have the same value (0) of "IL4"

FeaturePlot(scSeurat, feature = "IL6", split.by = "tissue_type") # Stronger in endo and mesem PMAIO

FeaturePlot(scSeurat, feature = "IL6-1", split.by = "tissue_type") # No signal
## Warning: All cells have the same value (0) of "IL6-1"

FeaturePlot(scSeurat, feature = "IL7", split.by = "tissue_type") # More in PMAIO cluster 4/5

FeaturePlot(scSeurat, feature = "IL10", split.by = "tissue_type") # More in PMAIO cluster 2

FeaturePlot(scSeurat, feature = "IL12A", split.by = "tissue_type") # No signal

FeaturePlot(scSeurat, feature = "TGFB1", split.by = "tissue_type") # Pretty even across all

Cluster 2 looks like it generally has more IL genes expressed. Isolate DE genes for this cluster.

clust1 <- sCVdata_list$res.0.4@DEvsRest$`2`[order(sCVdata_list$res.0.4@DEvsRest$`2`$Wstat,
                                                    decreasing = TRUE), ]
cat(rownames(clust1)[1:50])
## PTPRC TOX PRKCH INPP4B mikado.WCK01-AAH20201022-F8-rsc00190G443 ARID5B IKZF3 ARHGAP15 SATB1 DOCK10 TRPS1 RASA2 SAMSN1 NSD3 MBNL1 SH3KBP1 EMB SYTL3 ICOS mikado.WCK01-AAH20201022-F8-rsc00083G981 mikado.WCK01-AAH20201022-F8-rsc00275G164 THEMIS MSI2 SLAMF1 PITPNC1 NXF2;NXF2B;NXF5 SRSF7 IKZF1 RAPGEF6 STAT4 SLC38A1 RAP1A RABGAP1L ZC3HAV1 BCL2 TESPA1 APBB1IP GIMAP7 DHX8 DDX5 UBAC2 HUWE1 CD3E CTLA4 CD44 GIMAP1 RGS1 MTHFD1L RIPOR2 TANK

Visualise these genes in violin plots.

Interpretation: Now I see that IL2 and IL4 are expressed in more than just cluster 2, but it is nearly always in PMAIO. Cluster 8 is definitely where IL6 is expressed. IL7 is a bit more spread out across treatment types, but it’s definitely present in PMAIO in clusters 4, 5, and 11. Cluster 11 is plasma cells, but I’m honestly not sure of 4 and 5 - low signal hepatocytes? IL10 is in clusters 2 and 11 PMAIO. IL12A is a bit more in 4 and 5 PMAIO.

VlnPlot(scSeurat, feature = "IL2", split.by = "tissue_type") # Only in PMAIO cluster 2
## The default behaviour of split.by has changed.
## Separate violin plots are now plotted side-by-side.
## To restore the old behaviour of a single split violin,
## set split.plot = TRUE.
##       
## This message will be shown once per session.

VlnPlot(scSeurat, feature = "IL4", split.by = "tissue_type") # Only in PMAIO cluster 2

VlnPlot(scSeurat, feature = "IL6", split.by = "tissue_type") # Stronger in endo and mesem PMAIO

VlnPlot(scSeurat, feature = "IL6-1", split.by = "tissue_type") # No signal

VlnPlot(scSeurat, feature = "IL7", split.by = "tissue_type") # More in PMAIO cluster 4/5

VlnPlot(scSeurat, feature = "IL10", split.by = "tissue_type") # More in PMAIO cluster 2

VlnPlot(scSeurat, feature = "IL12A", split.by = "tissue_type") # A bit more in 1,4,5

VlnPlot(scSeurat, feature = "TGFB1", split.by = "tissue_type") # Pretty even across all